que invertida nos permite calcular los cuantiles (función cuantil):
\[Q(\text{p}) = F^{-1}(\text{p}).\]
Todo esto requiere integrar.
Para generar una intuición de estos conceptos, tomemos la posterior del modelo de número de incendios por verano. Aproximaremos las integrales discretizando los parámetros, una estrategia llamada integración numérica.
Integración numérica
Code
datos <-read.csv(here::here("R_curso_modelos", "datos", "barbera_data_fire_total_climate.csv"))# definimos parámetros de la previamu_a <-0; sigma_a <-1mu_b <-0; sigma_b <-0.1# Función que evalúa la densidad posterior no normalizada# ("un" de "unnormalized", porque siempre colonializado, nunca incolonializado)post_fire_un <-function(alpha, beta, log = T) { ## Log Verosimilitud (mismo código que en like_fire)# Calculamos la media, lambda: lambda <-exp(alpha + beta * datos$fwi)# usando lambda evaluamos la verosimilitud de cada observación, # en escala log: like_pointwise <-dpois(datos$fires, lambda, log = T)# la log-verosimilitud conjunta es la suma de las log-verosimilitudes por # observación: like <-sum(like_pointwise)## Log Previa lprior_a <-dnorm(alpha, mean = mu_a, sd = sigma_a, log = T) lprior_b <-dnorm(beta, mean = mu_b, sd = sigma_b, log = T) lprior <- lprior_a + lprior_b## Log Posterior lpost <- like + lpriorif (log) return(lpost)elsereturn(exp(lpost))}# Creamos grilla de valores de alpha y beta.side <-200# cantidad de valores por ladoalpha_seq <-seq(-0.75, 1.15, length.out = side)beta_seq <-seq(0.1, 0.22, length.out = side)grilla <-expand.grid(alpha = alpha_seq,beta = beta_seq)# Evaluamos la densidad posterior para cada valorsize <- side ^2# cantidad de valores en la grillafor (i in1:size) { grilla$post[i] <-post_fire_un(grilla$alpha[i], grilla$beta[i], log = F) }# Graficamos posterior no normalizadaggplot(grilla, aes(x = alpha, y = beta, fill = post)) +geom_tile() +scale_fill_viridis() +theme(legend.title =element_blank()) +xlab(expression(alpha)) +ylab(expression(beta)) +ggtitle("Posterior no normalizada") +nice_theme()
Calculamos la integral del denominador numéricamente, y obtenemos la posterior normalizada.
Code
size_a <-diff(alpha_seq)[1] size_b <-diff(beta_seq)[1]area <- size_a * size_b# área de cada celda en la grilla.grilla$post <-normalize_dens(grilla$post, area)# Graficamos posterior normalizadaggplot(grilla, aes(x = alpha, y = beta, fill = post)) +geom_tile() +scale_fill_viridis() +theme(legend.title =element_blank()) +xlab(expression(alpha)) +ylab(expression(beta)) +ggtitle("Posterior normalizada") +nice_theme()
Una vez que está normalizada, podemos marginalizar sumando (integrando.)
Code
# Para facilitar las cuentas, pondremos los valores de la posterior en una# matriz, análoga a la grilla (beta varía entre filas, alpha, entre columnas)post_mat <-matrix(grilla$post, side, side, byrow = T)alpha_tab <-data.frame(alpha = alpha_seq,post =NA)for (i in1:side) { alpha_tab$post[i] <-sum(post_mat[, i] * size_b)}plot(post ~ alpha, alpha_tab, type ="l", xlab =expression(alpha),ylab ="Densidad posterior")
alpha_tab$mass <- alpha_tab$post * size_a # los pesos, suman 1beta_tab$mass <- beta_tab$post * size_b# Esperanza de alphaea <-sum(alpha_tab$alpha * alpha_tab$mass)# Esperanza de betaeb <-sum(beta_tab$beta * beta_tab$mass)c(ea, eb)
[1] 0.1539278 0.1598357
Y cuantiles
Code
# obtenemos la función de distribución acumulada marginal para alphaalpha_tab$cdf <-cumsum(alpha_tab$mass)plot(cdf ~ alpha, alpha_tab, type ="l", xlab =expression(alpha),ylab ="Probabilidad acumulada")
Code
# creamos la función cuantil, es decir, la inversa de la acumulada.# ajusta una spline. Al tomar como x la probabilidad acumulada (cdf) y como # Y los valores de alpha, es la inversa de la acumulada (que es cdf = f(alpha))alpha_icdf <-splinefun(x = alpha_tab$cdf, y = alpha_tab$alpha, method ="monoH.FC") # esta función tiene como argumento la probabilidad acumulada (entre 0 y 1), # y devuelve el valor de alpha que acumula esa probabilidad (cuantil).probs <-c(0, 0.025, 0.5, 0.975, 1)quants <-alpha_icdf(probs)names(quants) <-paste(probs *100, "%")quants
Todo esto puede ser inviable si la posterior tiene \(> 3\) dimensiones (parámetros).
En algunos casos, la posterior tiene una solución analítica y además resulta en una distribución para la cual se conocen las marginales, con sus funciones de probabilidad acumulada y la inversa (función cuantil).
Pero eso es la excepción. Entonces, necesitamos formas más robustas de caracterizar las posteriores.
Enfoques para caracterizar distribuciones posteriores
Solución analítica con integrales manejables
Integración numérica
Aproximaciones determinísticas (Inferencia Variacional, INLA, Aproximación de Laplace)
Simulación
Simulación a partir de una grilla
Code
# (El mismo gráfico de antes)p1 <-ggplot(grilla, aes(x = alpha, y = beta, fill = post)) +geom_tile() +scale_fill_viridis() +theme(legend.title =element_blank()) +xlab(expression(alpha)) +ylab(expression(beta)) +ggtitle("Posterior normalizada") +ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +nice_theme()# Remuestreamos celdas de la grilla con reemplazo, donde la probabilidad es # proporcional a la densidad posteriornsim <-10000# número de muestrasrows_sim <-sample(1:size, size = nsim, replace = T, prob = grilla$post)# Tomamos las filas que salieron sorteadasdraws_table <- grilla[rows_sim, ]# Graficamos muestrasp2 <-ggplot(draws_table, aes(x = alpha, y = beta)) +geom_point(alpha =0.1) +ggtitle("Muestras de la posterior") +xlab(expression(alpha)) +ylab(expression(beta)) +ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +nice_theme()(p1 | p2)
A partir de las muestras es muy sencillo obtener las marginales y medidas de resumen:
Code
# Cada fila en draws_table es una muestra de la posterior conjunta (bivariada).# Si queremos ver la marginal de alpha o beta, simplemente miramos esa columna# en draws_table, ignorando los valores de la otra variable. La forma más rápida# de resumir una marginal es graficando su densidad empírica.par(mfrow =c(1, 2))plot(density(draws_table$alpha), xlab =expression(alpha), main =NA)plot(density(draws_table$beta), xlab =expression(beta), main =NA)
Code
par(mfrow =c(1, 1))# También podemos calcular medidas de resumen:summary(draws_table[, c("alpha", "beta")])
alpha beta
Min. :-0.72136 Min. :0.1012
1st Qu.:-0.02437 1st Qu.:0.1470
Median : 0.15704 Median :0.1597
Mean : 0.15440 Mean :0.1599
3rd Qu.: 0.33844 3rd Qu.:0.1724
Max. : 1.04497 Max. :0.2200
Code
# Para graficar ambas densidades en ggplot, alargamos:draws_long <-pivot_longer(draws_table[, c("alpha", "beta")], # sólo esas columnascols =1:2, names_to ="param", values_to ="value")ggplot(draws_long, aes(x = value, color = param, fill = param)) +geom_density(alpha =0.2) +facet_wrap(vars(param), scales ="free")
También podemos calcular la probabilidad posterior de cualquier tipo de afirmación sobre los parámetros o sobre funciones de los parámetros:
MCMC consiste en simular una cadena markoviana cuya distribución estacionaria es la distribución posterior.
Recorre el espacio de parámetros de forma estocástica, tomando muestras con una frecuencia proporcional a la probabilidad posterior.
Genera una secuencia de muestras de la posterior conjunta, generalmente correlacionadas.
Propiedades de las cadenas de Markov requeridas en MCMC
Irreducibilidad: desde cualquier punto debo poder llegar, en algún momento, a cualquier otro punto en el espacio de parámetros.
Aperiodicidad: no entra en ciclos infinitos.
Recurrencia positiva: en algún momento, podés volver a un lugar en dondes estuviste.
Bajo estas condiciones, luego de muchas iteraciones, las muestras tomadas corresponden a la distribución objetivo (la posterior), sin importar dónde comenzó la cadena.
Algoritmo de Metropolis (el MCMC más simple)
Defiinir el número de muestras a tomar (\(N\)), el valor inicial de los parámetros (\(\theta_1\)), y una distribución de propuestas (\(d\)) simétrica. Evaluar \(p(\theta_1 | y)\) (densidad posterior probablemente no normalizada).
Para \(i\) de \(2\) a \(N\), repetir:
Simular una muestra de la propuesta \(d\) centrada en \(\theta_{i-1}\). A esta muestra propuesta la llamaremos \(\theta^*\).
Aceptar la propuesta (\(\theta_i = \theta^*\)) con probabilidad \(\min(1, q)\). De lo contrario, definir \(\theta_i = \theta_{i-1}\).
Metropolis en R
Code
# Definimos la distribución de propuesta como una normal bivariada, y # aproximando la matriz de varianza-covarianza en base a la curvatura en el# modo. Para ello, optimizamos y obtenemos la hessiana.# optim necesita que la función a optimizar tome como argumento un vector # de parámetrospost_fire_opt <-function(x) { alpha <- x[1] beta <- x[2] lp <-post_fire_un(alpha, beta)return(lp)}opt <-optim(par =c(0 , 0), # vector inicial fn = post_fire_opt, # función a optimizarmethod ="BFGS", # método basado en gradientes (cool)control =list(fnscale =-1), # maximizar en vez de minimizarhessian =TRUE# que calcule la hessiana)# La matriz de varianza-covarianza, asumiento normalidad es la inversa de la # -hessianaV <-solve(-opt$hessian)# Para muestrear con una normal multivariada de propuesta, V puede escalarse# por un factor, que según la teoría, lo óptimo es 2.38 ^ 2 / d, con d siendo# la dimensión de la posterioroptimal_factor <-2.38^2/2# Creamos una función que genera una cadena markoviana para nuestra posterior.# start es un vector con el valor inicial de alpha y beta. n es el número de # muestras a tomar (incluye el start). factor es un factor que multiplica a # la matriz de varianza-covarianza de la propuesta; un buen valor es 2, pero es# interesante variarlo para ver cómo afecta la eficiencia.mcmc_fire <-function(start, n, factor = optimal_factor) {# Matriz para guardar las muestras draws <-matrix(NA, nrow = n, ncol =2) colnames(draws) <-c("alpha", "beta")# Vector para almacenar la log posterior, necesaria para comparar valores lp <-numeric(n)# Iniciamos la matriz de muestras y la lp draws[1, 1:2] <- start lp[1] <-post_fire_un(start[1], start[2], log = T)# Loop recursivo para hacer avanzar las coordenadasfor (i in2:n) {# Proponemos un nuevo valor (vector), simulando de la distribución de # propuesta, centrada en draws[i-1, ] y con matriza de vcov V. prop <-rmvn(1, mu = draws[i-1, ], V * factor)# Evaluamos la probabilidad posterior en la propuesta lp_prop <-post_fire_un(prop[1], prop[2], log = T)# Calculamos el cociente de densidad posterior propuesta / valor anterior.# En escala log, la resta equivale al cociente: lp_diff <- lp_prop - lp[i-1] p_q <-exp(lp_diff) # equivale a # p_q = exp(lp_prop) / exp(lp[i-1])# pero en log es más estable# Aceptamos si p_q > 1 (mayor densidad en la propuesta), y sorteamos que # la propuesta se acepte, con p_q si p_q es < 1. pkeep <-min(1, p_q) keep_prop <-ifelse(runif(1) < pkeep, T, F)if (keep_prop) { # si aceptamos la propuesta, la guardamos draws[i, ] <- prop lp[i] <- lp_prop # para comparar con el siguiente } else { # si no aceptamos, repetimos la muestra anterior draws[i, ] <- draws[i-1, ] lp[i] <- lp[i-1] # para comparar con el siguiente } }return(draws)}
Simulamos una cadena corta
Code
# Para graficar sobre la grilla, la extendemos un poco. Esto permite ver cómo la # cadena se mueve por zonas de baja densidad también, probablemente porque la # iniciamos en una zona de baja densidad.alpha_seq <-seq(-1, 1.5, length.out = side) # límites más amplios que antesbeta_seq <-seq(0, 0.3, length.out = side)size_a <-diff(alpha_seq)[1]size_b <-diff(beta_seq)[1]grilla <-expand.grid(alpha = alpha_seq,beta = beta_seq)for (i in1:size) { grilla$post[i] <-post_fire_un(grilla$alpha[i], grilla$beta[i], log = F)}
Code
# Corremos una cadena corta comenzando dentro de la grillarun1 <-mcmc_fire(start =c(-0.25, 0.15), n =50)draws_table <-as.data.frame(run1)# Graficamosggplot(grilla, aes(x = alpha, y = beta, fill = post)) +geom_tile() +scale_fill_viridis() +geom_point(data = draws_table, mapping =aes(x = alpha, y = beta),inherit.aes = F, size =2, alpha =0.3) +geom_path(data = draws_table, mapping =aes(x = alpha, y = beta),inherit.aes = F, alpha =0.5, linewidth =0.2) +theme(legend.position ="none") +xlab(expression(alpha)) +ylab(expression(beta)) +ylim(min(grilla$beta) - size_b, max(grilla$beta) + size_b) +xlim(min(grilla$alpha) - size_a, max(grilla$alpha) + size_a) +nice_theme()
# Corremos 3 cadenas de igual largo comenzando en puntos al azariter_warmup <-2000iter_sampling <-5000n <- iter_warmup + iter_samplingf <- optimal_factor *0.05c1 <-mcmc_fire(start =runif(2, -1, 1), n = n, factor = f)c2 <-mcmc_fire(start =runif(2, -1, 1), n = n, factor = f)c3 <-mcmc_fire(start =runif(2, -1, 1), n = n, factor = f)# Las Combinamos en un array 3D, ya que son 3 matricesarr <- abind::abind(list(c1, c2, c3), along =3)dimnames(arr) <-list(iteration =1:n,variable =c("alpha", "beta"),chain =1:3)# El paquete posterior requiere que las dimensiones sean iteration, chain, # variable, así que reordenamos:arr <-aperm(arr, c(1, 3, 2))# Converimos draws_array object, clase del paquete posteriordraws <-as_draws_array(arr)summarize_draws(draws)
# Quitamos el warmup, porque aún no había llegado al set típico (nos quedamos# con las últimas iter_sampling de cada cadena)iters_use <- (iter_warmup +1) : ndraws_crop <- draws[iters_use, , ]summarize_draws(draws_crop)
El \(\hat{R}\) (factor de reducción potencial de escala) verifica que las cadenas se hayan mezclado bien:
\(\hat{R} \to 1\) si hay convergencia.
Umbral recomendado: \(\hat{R} < 1.01\).
El número efectivo de muestras (ESS) indica a cuántas muestras independientes equivale tu muestra:
Considera la autocorrelación temporal en las cadenas.
Se recomienda ESS > 400 para estimaciones confiables (pero depende del contexto).
Definiciones:
Varianza entre cadenas (B de between): \[
B = \frac{N}{M-1} \sum_{m=1}^M \left( \bar{\theta}_m - \bar{\theta} \right)^2,
\] con \(N =\) número de muestras por cadena, \(M =\) número de cadenas, \(\bar{\theta}_m =\) media por cadena, y \(\bar{\theta} =\) media total.
Varianza intra-cadenas (W de within): \[
W = \frac{1}{M} \sum_{m=1}^M s_m^2,
\] siendo \(s^2_m\) la varianza de la cadena \(m\).
Varianza posterior marginal estimada: \[
\widehat{\text{var}}^+(\theta) = \frac{N-1}{N} W + \frac{1}{N} B
\]
Factor de reducción potencial de escala: \[
\hat{R} = \sqrt{ \frac{ \widehat{\text{var}}^+(\theta) }{ W } }
\]
Número Efectivo de Muestras: \[
\text{ESS} = \frac{MN}{\widehat{\tau}}
\]
donde \(\widehat{\tau}\) es una estimación de la autocorrelación temporal integrada en muchos lags.
Detalles en
Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. Rank-Normalization, Folding, and Localization: An Improved \(\hat{R}\) for Assessing Convergence of MCMC. Bayesian Analysis 16, Number 2, pp. 667–718.
Propuestas inteligentes: utilizando la curvatura de la densidad posterior (gradiente) logra dar pasos grandes manteniéndose en la región de alta probabilidad (set típico).
Para generar estas propuetas utiliza una analogía física: simulando el moviemiento de un cuerpo que se mueve en la superficie (N-dimensional) de la -log posterior siguiendo una dinámica Hamiltoniana.
Escala muy bien con la dimensionalidad (nro. de parámetros).
Como requiere calcular el gradiente, no se puede usar para estimar directamente parámetros discretos (hay que marginalizarlos).
Un programa de Stan define una log densidad que deseamos caracterizar.
En la mayor parte de los casos se tratará de una densidad posterior.
En este curso usaremos Stan para muestrear (HMC-NUTS), pero también ofrece otros métodos menos precisos pero más rápidos que valen la pena cuando tenemos sets de datos muy grandes o funciones de verosimilitud muy costosas.
El código
Stan divide el código en 6 secciones, pero rara vez usaremos todas
Code
functions {} data {}transformed data {}parameters {} transformed parameters {}model {} generated quantities {}
Ejemplo
Code
// (Esto se guarda en un archivo .stan)data {int N;array[N] int y; // Número de incendiosvector[N] x; // FWI}// Parámetros a muestrearparameters {real alpha;real beta;}// Calculamos las cantidades derivadastransformed parameters {vector[N] lambda = exp(alpha + beta * x);}// Acá definimos la log densidad posterior (o la que sea)model {// Densidad previa alpha ~ normal(0, 1); beta ~ normal(0, 0.1);// Verosimilitud y ~ poisson(lambda);}
Formas de definir la log densidad posterior
De forma expresiva, con sampling statements:
Code
model {// Densidad previa alpha ~ normal(0, 1); beta ~ normal(0, 0.1);// Verosimilitud y ~ poisson(lambda);/* Esto equivale a decir "sumá a la log densidad posterior la log densidad de alpha en una normal con media 0 y desvío 1, la log densidad de beta en una normal con media 0 y desvío 0.1, y la log verosimilitud (conjunta) de observar 'y' (es un vector) como realización de una Poisson cuyas medias son lambda (vector también). */}
De forma explícita, sumando términos a target:
Code
model{/* Stan guarda la log densidad objetivo en un objeto llamado 'target'. Cada vez que se evalúa la log posterior, target comienza valiendo 0. Esta parte de la función se encarga de sumarle los términos correspondientes. */// Sumar log densidad previatarget += normal_lpdf(alpha | 0, 1);target += normal_lpdf(beta | 0, 0.1);// Sumar log verosimilitudtarget += poisson_lpmf(y | lambda);// lpdf: Log Probability Density Function// lpmf: Log Probability Mass Function}
Y se pueden mezclar
Code
model{/* Podemos mezclar ambas formas, pero siempre hay que usar una sola por parámetro. (Si no, estará duplicada, lo cual es computacionalmente válido pero sería sospechoso que tenga sentido.) */// Vale hacer alpha ~ normal(0, 1);target += normal_lpdf(beta | 0, 0.1); y ~ poisson(lambda);/* Las funciones de densidad o verosimilitudes suelen tener constantes normalizadoras. Usando target +=, Stan las conserva, pero usando los sampling statements (~), los quita de la función, ahorrando un poco de cómputo. En algunos casos necesitamos conservar esas constantes para hacer una cuenta, y ahí conviene usar target +=. */}
Analogía con R
# Una función que evalúa la log posterior de nuestro modelo con un estilo # similar a Stanpost_fire_like_stan <-function(alpha, beta) { # Datos y <- datos$fires x <- datos$fwi# Parámetros# alpha y beta, pero en R no hay que declarar variables :) ... o :(, no sé# Parámetros transformados lambda <-exp(alpha + beta * x)# Modelo (definimos log densidad posterior) target <-0# iniciamos la log posterior en cero# Log densidad de las previas target <- target +dnorm(alpha, mean =0, sd =1, log = T) target <- target +dnorm(beta, mean =0, sd =0.1, log = T)# Es una analogía con el operador += de Stan# Log verosimilitud log_like_pointwise <-dpois(y, lambda, log = T) # dato por dato log_like_joint <-sum(log_like_pointwise) # conjunta target <- target + log_like_jointreturn(target)}
Interfaz entre R y Stan: cmdstanr
Code
# El código de Stan, como escribimos arriba, debe estar guardado en un archivo# .stan.# Compilamos el modelomodel <-cmdstan_model(here::here("R_curso_modelos", "modelos", "nfuegos.stan"))# Podemos revisar si lo escribimos bienmodel$check_syntax()# Creamos una lista nombrada para pasarle los datos. Los nombres de cada# elemento de la lista tienen que ser exactamente los que definimos en # la sección data {}stan_data <-list(N =nrow(datos),y = datos$fires,x = datos$fwi)# Y muestreamos la posteriorfit_mcmc <- model$sample(data = stan_data,chains =4, parallel_chains =4,iter_warmup =1000,iter_sampling =1000)
# Corremos los diagnósticos del muestreo. Si hay problemas, por defecto # aparecen warnings cuando sample() termina de correr, pero así podemos verlos # de nuevo.fit_mcmc$cmdstan_diagnose()
Processing csv files: /tmp/RtmpLFTFL1/nfuegos-202505070240-1-47a7ff.csv, /tmp/RtmpLFTFL1/nfuegos-202505070240-2-47a7ff.csv, /tmp/RtmpLFTFL1/nfuegos-202505070240-3-47a7ff.csv, /tmp/RtmpLFTFL1/nfuegos-202505070240-4-47a7ff.csv
Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.
Checking sampler transitions for divergences.
No divergent transitions found.
Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.
Effective sample size satisfactory.
Split R-hat values satisfactory all parameters.
Processing complete, no problems detected.
Code
# Resumimos las marginales usando el paquete 'posterior'summ <- fit_mcmc$summary()print(summ)
Algunas herramientas de visualización usando bayesplot
Code
# Seteamos un lindo color para bayesplotcolor_scheme_set("viridisC")# Extraemos las muestras y las formateamos como draws_array, una clase del# paquete 'posterior'draws <- fit_mcmc$draws(format ="draws_array")# Gráfico de traza, para visualizar convergenciamcmc_trace(draws, pars =c("alpha", "beta"))
Code
# Gráfico de densidad por cadenamcmc_dens_overlay(draws, pars =c("alpha", "beta"))
Code
# Sólo para ilustrar, reajustamos el modelo pero guardando las iteraciones# del warmup.fit_mcmc_all <- model$sample(data = stan_data,chains =4, parallel_chains =4,iter_warmup =1000,iter_sampling =1000,save_warmup =TRUE# por defecto es FALSE)
# Distribución predictiva posterior de y para FWI = 18:plot(density(ysim, from =0), xlab ="y", main =NA)
Code
# Visualizamos ambos a la vez:plot(density(lambda), xlab =expression(lambda, y), main =NA, xlim =c(0, 50))lines(density(ysim, from =0), xlab ="y", main =NA, col =4)
Code
# Y cuentas más arbitrarias aún. Ejemplo: ¿cuán probable es que el número de # incendios medio con un FWI = 20 sea el doble o mayor que cuando el FWI es 15?# Esto requiere calcular la distribución posterior de lambda para FWI = 15 y # FWI = 20, luego calcular el cociente y evaluar qué proporción de muestras # tienen un cociente >= 2.lambda_15 <-exp(d$alpha + d$beta *15)lambda_20 <-exp(d$alpha + d$beta *20)qlambdas <- lambda_20 / lambda_15plot(density(qlambdas), main =NA, xlab ="Cociente de lambdas")abline(v =2, lty =2)
Code
# Pr(qlambdas >= 2)sum(qlambdas >=2) / nsim
[1] 0.86225
Code
# También podemos reponder lo mismo sobre el número observado de incendios crudo,# no su media:y_15 <-rpois(nsim, lambda_15)y_20 <-rpois(nsim, lambda_20)qy <- y_20 / y_15plot(density(qy), main =NA, xlab ="Cociente de y")abline(v =2, lty =2)
Code
# Pr(qy >= 2)sum(qy >=2) / nsim
[1] 0.64375
Comparación entre previas y posteriores (marginales)
Code
par(mfrow =c(1, 2))plot(density(d$alpha), xlab =expression(alpha), main =NA, xlim =c(-3, 3))curve(dnorm(x, 0, 1), add = T, col =4)plot(density(d$beta), xlab =expression(beta), main =NA, xlim =c(-0.5, 0.5))curve(dnorm(x, 0, 0.1), add = T, col =4)